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1.  Introduction  and  Summary 

We  consider  a GI/G/l  queueing  system  in  which  customer  number  n , 

n=l,2,...  , arrives  at  time  x . We  let  x,  = 0 and  define 

n 1 


X 

n 


x - x 
n+1  n 


n 1 , 2 , • • • 


as  the  interarrival  time  between  the  nth  and  the  (n+l)st  customers.  We 
assume  that  X^,  X^,  ...  are  independent  identically  distributed  random 

variables  with  a cumulative  distribution  function  A(  ) , an  expectation 

2 

X , and  a variance  a . Let  M.  denote  the  number  of  the  customer  who 


initiates  the  kth  busy  period.  We  assume  that  the  distribution  of  the 
time  the  nth  customer  stays  in  service  is  a function  of  the  number  of 
customers  previously  served  in  that  same  busy  period;  namely,  denoting  by 


and  letting  j = n+1-2.  , then  the  nth  customer  is  the  jth  served  in  the 
busy  period  and  is  distributed  according  to  a distribution  G^ ( ) 


and  a variance 


/I" 


with  an  expectation 


00 
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Let  p,  = y,  / A . Then  if  there  exists  a k,  such  that  for  all 
k k 1 

k > k^  < 1 , the  queue  is  stable  and  the  stochastic  process  (M^; 

n=l,2,...}  is  positive  recurrent.  We  restrict  the  discussion  in  the 
present  paper  to  stable  queues. 

Denote  by  T the  time  spent  by  the  nth  customer  in  the  waiting  line, 
n 

Provided  that  the  distributions  of  X , n=l,2,...  , are  nonlattice,  then 

11  I 

T =>  T , where  =>  denotes  convergence  of  distributions.  Let  K denote 
n n 

the  number  of  customers  served  in  the  nth  busy  period  (K  = M ,,-M  ) . I 

n n+1  n 

Under  the  assumptions  stated  above  we  obtain  that  K^.K^,...  are  indepen- 

1 

dent  identically  distributed  random  variables.  In  Section  2 we  give  a 
numerical  procedure  for  the  determination  of  the  distributions  of  T 
and  K 

n 


Let  denote  the  length  of  the  nth  busy  period.  Clearly, 


B 

n 


Since  the  epochs  of  the  starting  of  busy  periods  are  regeneration  points, 
B^.B^...  are  independent  identically  distributed  random  variables. 


In  Section  3 we  present  a numerical  procedure  for  the  determination  of 
the  distribution  function  of  a busy  period. 


Several  papers  have  been  published  recently  on  numerical  analysis 
of  queueing  systems  (see  [1],  [3],  [4],  [5],  [6],  and  [7]).  Most  of 
these  algorithms  rely  on  the  memoryless  property  of  the  exponential  random 
variable.  In  the  current  algorithm  none  of  the  distributions  of  S^  or 

has  to  be  of  the  phase  type.  Furthermore,  both  distributions  can  de- 
pend on  the  number  of  customers  that  were  served  before  the  nth  customer 
in  the  same  busy  period.  The  disadvantage  of  the  current  algorithm  is 
that  one  must  execute  several  recursive  numerical  integrations.  Hence, 
when  traffic  is  heavy  the  computing  time  can  be  rather  long  and  the  results 
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may  contain  substantial  numerical  error.  Section  4 includes  a more  de- 
tailed discussion  of  the  results. 


2.  The  Distribution  of  the  Time  Spent 
in  the  Waiting  Line 


In  the  previous  section  we  denoted  by  T the  time  the  nth  custom- 

n 

er  spends  in  the  system;  we  also  noted  that  T — > T . In  this  section  we 

n 

determine  the  distribution  F of  the  random  variable  T . 


F(t)  = lim  P [T  > t]  , 
n-*»  n 


(2.1) 


and  denote  by  the  customer  who  initiated  the  busy  period  in  which 

the  nth  customer  is  served;  clearly, 

L = Max  {M.:  < n}  . 

n J 3 = 

To  calculate  P[T  > t]  we  use  the  fact  that  the  epochs  of  the  arrival  of 
n 


customers  to  an 


empty  system  (x^  , x^  , ...)  are  regeneration  points, 


hence  it  is  convenient  to  apply  the  following  conditioning, 


P [T  >t] 
n 


l P[T  >t  | L = n-j+1]  P[L  = n-j+1] 
j n n in 


We  first  notice  that  for  j=l  , 


P[T  >t  I L =n]  * 0 , for  all  t > 0 . 
n 1 n = 


When  j > 1 , 


P [T  >t  L 
n i 


(2.2) 


(2.3) 


. p[Y(Sn.j+1-Xn.1+1)>t  I Lo_j+1-n-j+l, 

'VjH'Vjtl1  * Jj)  Vi+1>>0] 

(2 

fj;1  j-l  1 

- p|  l (S  -X  )>t  I (S  -X  )>0,  ....  I (S  -X  )>o| 
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To  simplify  our  notation  we  let  H (t)  = P[S^-X^>t]  and 


H.(t)  = P 


r 2 j-2  j-i  ] 

(S  -X  )>0,  l (S  -X  )>0,  ....  I (S  -X  )>0,  l (S  -X  )>t  , 

L i=l  1 1 i=l  i=l  J 


so  that  expression  (2.4)  yields 

P [T  >t  I L =n- j+1] 
n 1 n J 


, t>0  , 


H (t)/H  (0)  . 


(2.5) 


Note  that  the  RHS  of  (2.5)  is  independent  on  n . As  for  the  calculation 
of  the  probability  P[L  =n-j+l]  , we  make  use  of  the  fact  that  (n-L^) 

is  the  backwards  recurrent  time  of  the  discrete  process  , j=l,2,... 

at  epoch  n . Hence,  denoting 


E[Kn]  = B , n=l , 2, . . . , 

and 

P [K  > k]  = P(k)  , 
n — 

we  obtain 

lim  P[L  =n-j+l]  = P(j)/B,  , (2.6) 

n-  n 


where  K (as  defined  in  the  introduction)  is  the  number  of  customers 
n 

served  in  the  nth  busy  period.  We  now  realize  that  the  numerator  of 
(2.6)  equals  the  denominator  of  (2.5);  thus,  combining  (2.1),  (2.2),  (2.5), 
and  (2.6),  we  get 

00 

F(t)  = j l H (t)  , t > 0 , (2.7) 

0 j=2  J 

where  B can  be  obtained  by  using 

OO 

B = 1 + l H (0)  . (2.8) 

j=2  J 

We  will  now  obtain  the  recursive  formulae  for  the  evaluation  of 
Hj (t)  . Let 
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T 
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U.(x)  = / G (u+x)dF(u)  , -«>  < x < ® , (2.9) 

J max{-x,0} 


then  the  calculation  of  ^(t)  is  straightforward: 

H2(t)  = P[S1-A1>t]  = 1 - Ux(t)  , t > 0 


(2.10) 


The  calculation  of  (t)  , j=3,4,...  , is  carried  out  recursively  as 
follows: 

°°  f 2 j-1 

H .+1(t)  = / P (S^)^,  l (S.-X  )>0,  ....  I (S  -X  )>0, 

J -°°  L i=l  i=l 

> (S.-X.)>t  I S.-X  =u  dP[S.-X.  < u] 

i=l  1 1 3 j J J 3 = 

00  r j-i  j-i  1 

= / P (S  -X  )>0,  ....  I (S.-X  )>0,  l (S.-X  )>t-u  dU  (u) 

-»  L i=l  1 i=l  J J 


/ H . (t-u)dU. (u)  + H. (0) (1-U  (t))  , j=3,4, . . . 

_oo  J 3 3 3 


(2.11) 


An  examination  of  expression  (2.31)  reveals  that  it  is  highly  unlikely 
that  tractable  analytic  expressions  for  H^(t)  , j=2,3,...  , can  be  ob- 
tained. We  are  compelled,  therefore,  to  approximate  F(t)  by 


N / N \ 

,(t)  = ^ Hj  (t)  / ll  + ^ Hj  (0)1  , 


(2.12) 


while  the  Hj (t)  , j=2,3,...  are  calculated  numerically. 

To  assess  the  magnitude  of  the  numerical  error  that  results  from 

summing  up  only  N elements  in  (2.12),  one  should  first  assess  the  value 
-00 

of  4j_N+1  Hj (t)  . This  can  often  be  achieved  by  using  the  fact  that 


Hj (t)  < P 


E <srt’  ’ c] 


(2.13) 


T-372 


The  RHS  of  (2.13)  can  sometimes  be  approximated  using  the  central  limit 
theorem.  If  the  c.d.f. 's  , k=l,2,...  , satisfy  < G^  , where  G^ 

2 

is  a c.d.f.  having  an  expectation  y^  (y^<A)  and  a variance  , then 


for  j sufficiently  large  we  may  write 

= 1 - <t> 


>1  1 
5 l (S  "X  )>t 
_i=l 


ft  - (j-l>(U0-A)> 


(vo + °2), 


where  <J>  denotes  the  standard  normal  c.d.f.  If  (t  - ( j-1) (y^-A))  / 
j-1) | Vq  + <?2 j > 1 , then  after  some  algebraic  manipulations  we  obtain 

2> 


't  - (j-i)(yn-A)  x , 

1 - 4>l  — \ < exp|-t 

/2tt 


X - u 


0 


vV(J'1)(vo  + °2), 


v2  + a2 ' 
0 


exp  - 


(j-l)(A-y0)‘ 
+ a2' 


V + °2)  /' 


Now  we  apply  the  facts  that  H.(0)  > H 1 (0)  and  H.(0)  > H.(t)  , and 

v J = 3+1  j = 3 

denoting 

2> 


N^_  = min 


1 / ^ ^0 

j : «N(0)  > exp|-t  — 2 — 2 • exPl“ 

V^TT 


j=N+l , N+2, . • • > , 


we  get 


V0+O 


(j-l)(X-y0)' 


l H (t)  < C (t) 

j=N+l  3 


(2.14) 


where 


CN(t)  = <Nt-N)HN(0)  + 


exp) 

-a-yQ)(t  - 

O-M0>»t/2)  / (vj«2)) 

*^2tt  | 

|l  - expj 

[-(A-Uo)2  / 2|v2«2|) 

) 

Combining  (2.7),  (2.8),  (2.12),  and  (2.14)  yields 
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N 

l 


J-2 


I 

Jbl 


H . ( t) 
3 


H.(t)  + CN(t) 


VO 

F(t) 


1 + 


y°> 


N 

1 + l H (0) 


(2.16) 


The  cases  in  which  use  of  the  central  limit  theorem  cannot  be  justified 
should  be  handled  on  an  individual  basis.  Here  it  may  sometimes  be  very 
difficult  to  obtain  the  required  bounds. 


Numerical  examples 

A computer  program  was  written  for  the  calculation  of  Fg(t)  . 

The  recursive  integrations  (2.11)  are  carried  out  via  the  discretization 
of  X and  S . The  ranges  for  which  the  densities  are  substantial  are 
divided  into  intervals  of  length  A . It  is  assumed  that  the  discretized 
random  variables  take  the  value  of  the  middle  of  an  interval  with  a prob- 
ability equal  to  the  concentration  of  their  densities  in  the  interval. 


Example  1 : Interarrival  times  gamma  distributed. 

0 < x . 


dA(x) 

dx 


7.010'5  9.5  -7. Ox 

x e 


r(10.5) 


The  service  times  are  normally  distributed  ~ = 0 . 95  |l  + ^ 3 j , 

0-0. 15 J . The  function  F^(t)  for  N=7  and  A=0.025  (computing  time 

40  seconds  on  the  IBM  370/148)  is  given  in  Figure  1.  The  RHS  and  LHS  of 
Equation  (2.16)  were  also  calcula 

the  results  are  given  in  Table  1. 


Equation  (2.16)  were  also  calculated  taking  y^  = 1.0  and  = 0.15  ; 


TABLE  1 


THE 

RHS  AND 

LHS  OF 

(2.16) 

FOR 

THE  DATA 

OF  EXAMPLE  1 

n 

0.0 

0.1 

0.2 

0.3 

0.4 

0.5 

0.6 

0.7 

0.8 

RHS 

1.00 

1.00 

1.00 

1.00 

1.00 

1.01 

1.01 

1.01 

1.01 

LHS 

0.98 

0.98 

0.97 

0.96 

0.93 

0.90 

0.83 

0.74 

0.62 
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Example  2:  Interarrival  times  uniformly  distributed  on  (0,1]. 
Service  times  have  the  following  distribution: 


dGk(x) 


1 

2 


1 

2 


(2k-l) ! k-1 

; 2 x 

(k-1)  ! 


for  k>4  we  have  G,  = G,  . 

k 4 


(l-x)^  ^ dx  , 
The  function 


x=0 

0<x<l 

Fq  (t)  for 


N=15 


k-1, 2, 3,4  ; 


and  A=0.025 


(computing  time  52  seconds  on  the  IBM  370/148)  is  given  in  Figure  1.  The 
bounds  (2.16)  were  calculated  here  taking  Gq  = G]  ’ anc*  t*le  resu-^ts  are 


presented  in  Table  2. 


TABLE  2 


THE  RHS 

AND  LHS 

OF  (2. 

16)  FOR 

THE 

DATA  OF 

EXAMPLE  2 

0.0 

0.2 

0.4 

0.6 

0.8 

1.0 

BBI 

RHS 

1.03 

1.03 

1.03 

1.04 

1.04 

1.04 

1.05  1.06 

LHS 

0.93 

0.91 

0.88 

0.83 

0.74 

0.63 

0.43  0.32 

3.  The  Distribution  of  the  Length 
of  a Busy  Period 

In  the  introduction  we  denoted  by  B the  length  of  the  nth  busy 

n 

period  and  noticed  that  ®i’B2»*.‘  are  independent  identically  distributed 

random  variables.  In  this  section  we  calculate  the  cumulative  distribution 

function  of  B . Define 
n 

R(t)  = P[B  >t]  , t > 0 . 
n 

We  calculate  R(t)  using  the  following  formula: 


R(t)  = y P[B  >t,  K =k]  . 
k=l  n n 


(3.1) 


T-372 


Denote 


and 


Q1(t,x)  = P[S  >t,  S1-A1>x] 


Qn(t,x)  = P 


[n  n- 1 n 1 

l S >t,  S -A  >0,  I (S  -A  )>0,  l (S  -A. )>x  , 

i=l  1 i=l  11  i=l  1 1 J 


n=2, 3, ...  . 

The  elements  in  (3.1)  are  calculated  as  follows.  For  k=l  we  obtain 

P[Bn>,:’  V1]  = f I1“F1^S>)  dGiCs)  = 1 - G1(t)  - Q1(t,0)  ; (3. 

for  k>l  we  get 

pk>t‘  Kn=kl  = p l si>t:’  SrV° Tcs.-A.)^,  I (S.-A  )<o| 

J Li=l  i=l  i=l  1 1 J 

r k k-i  ] 

l ^>t:»  S -A  >0,  ....  I (S  -A  )>0 

U=1  1 i=l  1 1 J 

(ii Si>t'  Si"v0,  •••’  il1<si_Ai)>°] 

” [k  k-1  "I 

-0  P|^l  Si>t’  W0’  ^ (Si-A1)>0  | Sk=s  dGk(s) 


2) 


= P 


(3.3) 


- Qk(t,0) 


= fQ  Qk-l(t"S’  0)dGk(s)  + Qk_iC°,0> (i-GkCt))  - Qk(t,0) 
The  functions  Q (t,x)  , k _>  1 , can  be  calculated  recursively.  First, 


Q.(t,x)  = / F (s-x)dG  (s)  , 

1 max{t,x)  1 1 


(3.4) 


then  for  k > 1 we  have 


10  - 


T 
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oo  oo  f*  k k-1 

Qk(t,x)  = / / P I S^t,  S -A  >0,  ....  I (S.-A.)>0, 

y=0  s=0  Li=l  i=l 


i ‘W*  l Vs 


i=l 


* Vy] 


dGk(s)dF(y) 


oo  oo  j"k-l  k-1 

/ / p I S.>t-s,  S -A  >0,  ....  I (S.-A  )>max{0,  x-s-l-y} 

r=0  s=0  l_i=l  1 11  i=l  1 1 


dGk(s)dF(y) 


After  some  manipulations  we  obtain  for  x < t , 
t-x  x+y 

Q,,(t,x)  = / f Q^-^t-s,  x-s+y)dG,  (s)dF(y) 


y=0  s=0 


x+y 


+ f f Q,i(0,  x-s+y)dG  (s)dF(y) 
y=t+x  s=0 

t-x  t 

+ f f Qk_1(t-s,  0)dGk(s)dFk(y) 

y=0  s=x+y 


(3.5) 


+ Q 


1(0,0)  (l-Gk(t)F(t-x)|  + / (l  - Gk(x+y))  dFk(y) 

L y=t-x  J 


and  if  x > t then 


Qk(t,x)  = Qk(x,x)  . 


As  in  the  previous  section,  we  approximate  R(t)  by 

N 

P f B >t.  1 

n 


(3.6) 


Rn(t)  = l P[B  >t,  K.-k] 
k=l 


(3.7) 


To  assess  the  sum  of  the  remaining  elements  we  can  use  the  fact  that 
for  k > 1 , 


l K =k]  < l P[K  =k]  = H (0) 

k=N+l  n k=N+l  i 


- 11  - 


In  Figure  2 we  present  the  function  R^ft)  for  the  two  examples 

specified  in  Section  2.  The  numerical  integration  is  carried  out  here, 
as  in  the  examples  of  Section  2,  via  the  discretization  of  the  random 
variables.  This  program  requires  much  longer  computing  periods  because 
two-dimensional  integrations  have  to  be  carried  out  in  many  points.  The 
computation  of  Rg(t)  on  t*le  370/148  for  the  data  of  Example  1, 

taking  N=4  and  A=0.15  , required  six  minutes  (H^ (0)  = 0.011)  , and 

the  computations  for  Example  2's  data  (N=5  , A=0.1)  also  required  six 
minutes  (H^(0)  = 0.046)  . 

4.  Conclusions 

This  paper  presents  numerical  procedures  for  the  determination 
of  the  distributions  of  the  time  a customer  spends  in  the  waiting  line  and 
the  distribution  of  a busy  period.  The  algorithm  can  use  standard  distri- 
bution functions  as  well  as  empirical  distribution  functions,  and  it  can 
efficiently  handle  service  times  and  interarrival  times  that  are  discrete 
random  variables.  Although  our  computer  programs  are  quite  inefficient 
and  we  believe  that  the  computing  times  can  be  substantially  reduced,  it 
is  expected  that  the  required  CPU  times  will  generally  be  long,  especially 
in  heavy  traffic.  When  the  distributions  involved  are  continuous,  the 
smaller  A is,  the  more  accurate  are  the  results,  and  the  longer  are  thei 
computing  times.  The  "bounds"  given  in  (2.16)  are  only  approximate  bounds 
since  they  make  use  of  the  central  limit  theorem.  It  is  expected,  however 
that  they  give  a good  indication  of  the  error  resulting  from  replacing 
F(t)  and  R(t)  by  F^Ct)  and  Rq(0  ’ resPect^vely* 
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